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Peebles (1989) showed that in the gravitational instability picture galaxy orbits can be traced back in 
time from a knowledge of their current positions, via a variational principle. We modify this variational 
principle so that galaxy redshifts can be input instead of distances, thereby recovering the distances. As 
a test problem, we apply the new method to a Local Group model. We infer M = 4 to 8 x IO^^Mq 
depending on cosmology, implying that the dynamics of the outlying Local Group dwarves are consistent 
with the timing argument. Some algorithmic issues need to be addressed before the method can be 
applied to recover nonlinear evolution from large redshift surveys, but there are no more difficulties in 
principle. 

Subject headings: large-scale structure of universe - galaxies: distances and redshifts - Local Group 



1. INTRODUCTION 

Phase space is six-dimensional and therefore six num- 
bers for each particle will specify the dynamics completely. 
The standard approach is to define initial conditions as six 
numbers (positions and velocities in three dimensions) and 
integrate those forward in time. This is the usual V-body 
approach to the problem. 

Alternatively, it is not necessary to define those six num- 
bers at only one time. It is equally well possible to split 
them, such that three numbers will be known at an initial 
time, and three at a final time, where the former are de- 
rivcd from physical arguments about the inital atatc, and 
the latter are provided b}^ data on the current etato of the 



system. This is the boundary value approach. 

The usual boundary value for structure formation by 
hierarchical clustering at initial times is the gravitational 
instability requirement that initial peculiar velocities must 
vanish: therefore, it is possible to express the orbits as a 
sum of growing modes. This is what perturbation the- 
ory is designed to do, and linear perturbation theory and 
its extension, the Zel'dovich approximation are in wide 
use. However, as structure formation is non-linear on small 
scales, a non-linear formalism would be more useful. Non- 
linear perturbation theory exists (e.g. Nusser and Deke 
(19931, iGramann (1993a)| , |G ramann (1993b)|, p3uchert ancj 



Ehlers (1993)) but is very complicated and has conver- 
gence problems, i.e. the dynamics at early times can be 
fitted to a high accuracy only at the expense of a good fit 
at later times. 



A different method (cf Peebles (1989)) of addressing the 
boundary value problem is to use the variational princi- 
ple. The basic method is to start with a parameterisation 
of the orbits which satisfies the boundary conditions and 
then adjust the parameters until a statio nary point of the 
action is found. (A variant, suggested by Giavalisco ct al 



1993 and implemented by Susperregi and Binney (1994 



parameterises the density and velocity fields rather than 
the orbits.) This may seem like perturbation theory be- 
cause it is a method which attempts to get better and 
better approximations of galaxy orbits but there are im- 
portant differences. The main one is that perturbation 



theory attempts to fit early times even at the expense of 
accuracy at later times, whereas the variational principle 
spreads out the errors more uniformly over all times. The 
variational principle is also algorithmically more straight- 
forward to do to higher orders. 

The main disadvantage of Peebles' original method is 
that it requires the input of distances to recover red- 
shifts. Redshifts, however, are easy to measure, whereas 
distances are not. We therefore change to recovering dis- 
tances from redshifts. The recent nearly all-sky redshift 
surveys QDOT and PSCz provide a s trong motivation for 
variational m etho ds in redshift space. | Shaya, Peebles, aiu^ 



Tully (1993) and 3haya, Peebles, and TuUy (1995) have 



pointed out that this could be achieved by treating the 
distance boundary condition as a parameter and fitting 
that until the recovered redshifts agree with the measured 
ones. Another approach is to modify the variational prin- 
ciple until the boundary conditions are of the desired fo rm. 
We take such an approach and so does Whiting (1998)| but 
the details in his treatment differ from ours. The attrac- 
tive feature of this approach is that it only requires small 
modifications of Peebles' elegant original method. 

In this paper we develop a redshift space variational 
method, apply it to a small system — the Local Group 
— , and point out the algorithmic problems that need to 
be solved to extend to larger systems. 

Various aspects of the Local Group dynamics have been 
stud ied by [Peebles (19891, [Peebles (1990)| , [Peebles (1991) , 
and Dunn and Laflamme (1993)|. The new feature of our 



analysis is that we attempt to constrain the masses of the 
Milky Way, M31, and an underlying distribution of un- 
clustered matter by using a likelihood approach. 

Two points need to be made about these variational 
methods in general: Firstly, the true solutions of the vari- 
ational action need not be a minimum eve n though the co l- 
loquial use of 'least action' has stuck (cf Peebles (1990) ). 
Secondly, the equations solved are the same as for A^-body 
integration, only the approximations used to solve them 
are different. In neither of the two approaches do parti- 
cles have to be galaxies - they may well be s amplers of 
some underlying distribution function and as Branchin: 
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and Carlberg (1994) and Dunn and Laflamme (1995) have 
pointed out, applying the variational method using galax- 
ies only, neglecting biasing and unclustcred matter, leads 
to wrong results. 

2. FORMULATION 
Orbits of galaxies follow equations of motion of the form 



-(a'i.) = - — 



ax 



where x are some comoving coordinates, and a{t) is the 
scale factor of the universe. These equations are deriv- 
able from a stationary action, 5S — 0, subject to some 
boundary conditions. In this work, we will distinguish be- 
tween two different actions and different sets of boundary 
conditions leading to the same equations of motion. 

2.1. Two Kinds oj Boundary Condition 

For simplicity, let x be one-dimensional in this subsec- 
tion. Then the two cases are: 

• the 'real space' case, where the boundary conditions 
are a^x = at t = (i.e., initial peculiar velocities 
zero) and 5x — Q aX t — \ (i.e., present positions 
fixed) with the action given by 



S = 



Ldt 



where L = -^a x 



- H^) (2) 



This is the action and the boundary conditions first 



proposed by Peebles (1989) 



• the 'redshift space' case, where the initial boundary 
condition remains the same and the final one be- 
comes aSx + aSx = at t = 1 (i.e. present redshifts 
are fixed). The new boundary conditions will still 
lead to the same equations of motion if we change 
the action to 



L-G]dt , where G 



a^xx - 



^aax^ 



(3) 



and L remains the same. Whiting (1998) describes 
more general transformations of this type. 

The standard method for finding an orbit that will keep 
the action stationary and therefore is a solution to the 
equations of motion consists of formulating a parametric 
expression for the coordinates x which satisfies the bound- 
ary conditions and adjusting the parameters to make the 
action stationary. For the real space case, we choose the 
expression 

N 

X = Y,CnUt), (4) 

n=l 

where C„ are a series of N coefficients and fn{t) is a tem- 
poral basis function. For the redshift space case. 



cz 

W) 



N 

E 



C„5„(<) 



(5) 



To ensure that the boundary conditions are fulfilled, the 
basis functions have to satisfy 



/n(l) = 



hit) 



a(l)' 



9n>l{t) = fnit) (6) 



A good choice for /„ is (1 — D{t))'^, where D{t) is the 
growth factor from linear theory; this was proposed by |Gi-| 



avalisco et al. 1993 , In this case, for iV = 1 the variational 
method reduces to linear theory. But all that is essential 
is that linear theory should be followed as i ^ 0. So, other 
choices such as /„ = (1 — t^/'^)" are also possible. 

2.2. Equations for Stationary Action 

The real three-dimensional problem is a combination 
of the two types of boundary conditions. Because of the 
way in which we choose coordinate systems, the real space 
treatment can be used for the first two dimensions {x, y) 
of each orbit and only the third dimension (z) will have 
to be treated in redshift space. We will refer to axes by 
their dimension rather than by a x, y, or z in order to 
prevent confusion between the third dimension and the 
symbol used for redshifts. 

One coordinate system is assigned to each object and 
all coordinate systems share a common origin. This ori- 
gin is the point from which a comoving observer would 
measure the objects' redshifts. It is important to note 
that, although the objects move, their coordinate frames 
do not. By pointing the 3- axis of each frame towards the 
current galactic coordinates l,b of the object associated 
with it, we ensure that the radial velocity is along one axis 
only. Therefore, the redshift space treatment will have 
to be applied only to that axis. With this orientation, 
the 3-coordinate of each object is czi (czi = dxi^ + axi^), 
where Zi is the 'comoving' redshift, i.e. the redshift mea- 
sured by an observer who is at rest with respect to the 
microwave background. The 1,2-coordinates at t = 1 are 
— 0, Xirj — 0, so the current positions are completely 
specified and real space treatment can be applied. 

In these coordinates, the parametric expansion of the 
position X of an object i becomes 



N 



N 



Jnit) 



1(1) ^ 

^ ^ n— 1 

and the full action is Jq{L — G), where L and G are 
L ^ ^ m.iax^--^aa 2_^mi-K^-\ — 2^ 



and 



G ^^m, (a^Xi^^^x.^^^ + ^aaxf^^ 



(7) 
(8) 

'i'tidal 

(9) 
(10) 



The first term in equation |^ describes the total kinetic 
energy, the second term is an acceleration caused by the 
fact that the coordinate system is expanding at a rate a{t), 
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and the third term describes the gravitational interaction 
between members of the group of objects. $tidai repre- 
sents the tidal potential caused by the influence of objects 
external to the region considered. Note that introducing 
a homogeneous mass distribution into the system is very 
easy since it will only rescale the second term. 

Inserting the equations ^ and |l^ into the action 5* and 
taking the gradient with respect to the Ci leads to 



a frnfndt 



^i(l,2): 



Jo a C'*-'i(l,2) 



(11) 



and 



/i(i)V 
d(l) 



aagndt—— 
a(l) 



a gragndt 



g-n \ ^ 
a ^ 



m,mi\7i,,, |x,; - xJ ^ 



tidal 



da 



(12) 



H3) 



To find the orbits which keep the action stationary, we 
now have to solve equations |ll| and |l^ for the spatial coef- 
ficients Ci^n- The standard method is to assume some ini- 
tial coefficients, calculate the right hand sides, then solve 
the left hand sides for the coefficients and use the results 
to recalculate the right hand sides. Note that the left hand 
side integrals in equations |ll| and |l^ have to be done only 
once. 

We have experimented with various procedures to 
help convergence, e.g. adding a stabilising term oc 
[/ da frnfndt] Ci^m On both sides of the equation or scaling 
down some terms for early iteratio ns and only slowly in- 
creasin g them to their full val ue (cf Susperregi and Binney 
(1994|l, ]Giavalisco et al. 1993|) . 



2.3. Units 

It is convenient to define model units for time, mass, 
and length (tick, marble, and lap) which we keep scaleable 
with the age of the universe Tq, since Tq is not known. The 
units are given by 



1 tick = Tn = 10 K Ga 



(13) 



3. APPLICATION TO THE LOCAL GROUP 

In this section we will apply the above formalism to 
a small group of galaxies, namely the Milky Way-M31 
system and outlying Local Group dwarves, with external 
forces approximated by dipole and quadrupole tidal forces 
growing according to linear theory. Since distances to ob- 
jects within this group are fairly well known, we can con- 
strain the mass distribution by comparing the distances 
predicted by our method to the observed distances for a 
range of different mass distributions. We also test the ef- 
fects of assuming different values for fio . 



3.1. Formulation for the Local Group 

The formalism is outlined in section ^, but for the Local 
Group there are two special considerations: 

Firstly, the Milky Way has a special role. As explained 
above, each galaxy has its own coordinate system, and we 
take the current position of the Milky Way as the com- 
mon origin. This origin is at rest with respect to the CMB 
and for t ^ 1 the Milky Way moves away from it. As 
the Milky Way's current position is fixed the real space 
treatment can be applied in all three dimensions. 

Secondly, cosmology only influences the Local Group via 
tidal forces. Therefore, we can keep things simple by in- 
troducing the following complication: The integrals on the 
left hand sides of equations |ll| and |l^ are particularly 
simple to solve for a{t) = t^^^, and we exploit this fact by 
setting a{t) = t^/^ regardless of cosmologies. This a{t) is 
only the scale factor of a convenient frame that we have 
chosen for the Local Group and which need not necessar- 
ily bear any relation to the real physical scale factor of the 
universe aphys(i) except at the boundaries 



a{t ^ 0) = t^/^ cx aphys(i 0) 
a(t = l) = 1 = aphys{t = 1) 



(16) 



For interactions between members of the group, the value 
of a{t) does not matter. Only when calculating the in- 
fluence of the tidal forces do we have to consider the 
relation between the two sets of comoving coordinates 



Xphys 



-X. The tidal potential then is 



tidal 



m 

^phys 



[d ■ Xp];^ys + 2"^ptiys ' Q ' -^phys] (l'^) 



and K is the age of the universe in units of 10 Ga. We 
do not want the velocities to be scaleable, so we require 
100 km/s = 1 lap/tick which leads to 



1 lap = 1.023 K Mpc 



(14) 



The assumption of a gravitational constant equal to unity 
as in the equations above leads to 



1 marble = 2.38 ■ 10^^ k M, 







(15) 



We therefore have a one-parameter family of models for 
different ages of the universe. 



where D{t) is the growth rate of the density fluctuations 
d is the dipole and Q is the quadrupole st rength 



The quadrupole poten tial is taken from Raychaudhury 
and Lynden-Bell, (1989)| . The dipole is inferred indirectly 
from the Milky Way's known velocity with respect to the 
microwave background, in the following way. As the equa- 
tions (11) and (12) are being solved for the orbits of all 
galaxies involved, we constantly monitor the development 
of the predicted velocity of the Milky Way. This velocity 
depends on the strength of the dipole. During the itera- 
tions, we fit that dipole strength such that the predicted 
vel ocity of the Milky Way agrees with the one measured 
by |Kogut ct al. 1993| . 
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3.2. Finding Solutions in the Local Group 

We model the Local Group as an ensemble of 13 galax- 
ies with the mass of the system distributed between the 
Milky Way, M31,and some local unclustered matter psm- 
It is convenient to express this Psm in units of the critical 
density, but it is important to note that it is only local and 
does not change the cosmology. The cosmological models 
have A = and various fio- The ratio of masses of Milky 
Way and M31 is fixed at 2:3 and all dwarf galaxies are 
treated as test particles. 

We select the dwarf galaxies by cho osing all objects from 



the li pt of Local Group mcmbcra (Hodge ct al (1003) ^ 



whichi are more than 500 kpc away from both M31 and tho 



Milky Way. This high distance is a necessary restriction 
because the orbits of nearby galaxies tend to be too dom- 
inated by the internal dynamics of the Milky Way-M31 
system and are therefore difficult to reconstruct. Nearby 
galaxies might also be part of the system's halo. The dwarf 
galaxies chosen, on the other hand, are generally not too 
sensitive to the mass ratio between M31 and the Milky 
Way because of their r elatively great distan ce to that sys- 
tem. Distances (from [Hodge et al (1993) ) and redshifts 
(from the Nasa/IPAC Extragalactic Database) of all galax- 
ies are listed in table 1. Note that these heliocentric red- 
shifts are converted to CMB-centric redshifts to get the 
appropriate boundary conditions. 



Galaxy 


/ 


b 


cz 


^^obs 


f^modcl 




deg 


deg 


km 

s 


kpc 


(kpc/s) 


M31 


121.2 


-21.6 


-300 


725 


787 


IC1613 


129.8 


-60.6 


-234 


765 


498 


WLM 


75.9 


-73.6 


-116 


940 


1224 


Sextans A 


246.2 


39.9 


324 


1300 


1878 


N3109 


262.1 


23.1 


403 


1260 


1977 


ICIO 


119.0 


-3.3 


-344 


1250 


1227 


Pegasus 


94.8 


-43.5 


-183 


1800 


1911 


Sextans B 


233.2 


43.8 


301 


1300 


2023 


SagDIG 


21.1 


-16.3 


-77 


1150 


1266 


LGS 3 


126.8 


-40.9 


-277 


760 


1220 


EGB0427-h63 


144.7 


10.5 


-99 


800 


2069 


N6822 


25.3 


-18.4 


-57 


540 


1554 



Table 1: Local Group Members; predicted distances are from 
the solution pictured in figure 1 

In order to increase our chances of finding the right so- 
lutions among the many possible ones, we start the code 
with an initial guess by choosing the Ci^„ = except for 



C'io.i — di 



CZj 



(18) 



This produces an initial set of Xi{t) which fit the red- 
shifts and measured distances (if 1 lap = 1 Mpc) but 
which are not solutions of the equations of motion. We 
therefore add an extra term to equations ^ and ^ to 
change them into something that the are solutions of. 
This term is then gradually removed from the equations 
during iteration. 

Our current method always produces convergence of the 
coefficients, but does not do so in a very efficient way. In 
most cases, we need several hundred iterations to achieve 
convergence, where the main problem seems to be that 
some of the converging Ci^„ are caught in a cycle of two 



values. We have tried to remedy this problem in a very 
crude way and were successful but only at a severe cost 
in iterations; this is one of the problems that will have to 
be solved in a more general way, before the code can be 
apphed to large datasets. 

We check our reconstructed by taking their values at 
t = 0.01 as the initial conditions of an A^-body integration. 

Figure 1 shows the variational orbits and A^-body orbits 
for a system of a total mass of 3.10 marbles. For the mass 
distribution of table 1, the dipole acceleration inferred is 
0.85 laps/tick^ in the direction I = 274° and b = 29° (the 
direction is indi stinguishable from Yahil, Tammann, ancj 



Sandage (1977)) 




Figure 1. Variational and A-body orbits for Q = 0.99, 
mtot = 3.10 marbles (no unclustered matter); squares indicate 
positions a,t t — 1; the scale is in comoving laps and the coordi- 
nate directions are (cos 6cos Z, cos 6sin Z, sin 6). We have labelled 
as many of the galaxies as was possible without overcrowding 
the figure. 

The last column of table 1 lists the predicted distances 
for this particular model. They are generally similar to the 
measured distances, which is not a trivial result, since as 
far as the mathematical formalism is concerned, the pre- 
dicted distances are not even constrained to be positive. 
In fact, when trying to reconstruct the orbits of galaxies 
which are too close to M31 or the Milky Way, the code does 
produce negative distances. A typical example of this is 
Leo I (cf ^aritsky ct al (1989) ): the code is given a pos- 
itive redshift but only finds approaching solutions. The 
only way to reconcile these two requirements is to put the 
dwarf galaxy on the other side of the Milky Way - which is 
what produces the negative distance. The reason for this 
is that nearby galaxies are too dominated by halo dynam- 
ics and we have not attempted to model the halo in any 
way. Hence our decision to exclude all the nearer dwarf 
galaxies. 

3.3. Analysis of Local Group Solutions 

To quantitatively analyse our results, we assign a like- 
lihood to each set of solutions and therefore to each mass 
distribution. The following is very crude, but we decided 
to use it because it is well defined and uses plausible 
reasoning. The likelihood is calculated in the standard 
Bayesian fashion by assessing the probability of a set of 
solutions given the set of measured distances. We have to 
allow for the uncertainties associated with the measured 
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distances, and the fact that some of our solutions might 
simply be wrong (we call this the 'outlier probability'). As 
we have no particular reason to believe that the uncertain- 
ties in dohs are Gaussian, we use (for all galaxies except 
M31) the hatbox function 



hb(a;, a, b) 



ib- 




a)- 



\i a < X <h 
otherwise 



If the outlier probability is a and the uncertainty in dobs 
is /3 then 

prob((iobs I c'modei) = ahb (dobs, min((iobs), max(dobs)) 
+(1 - a)hb ((ij,bs, (1 - /3)<bs, (1 + /3)<bs) (19) 

For M31 we take prob((i*|_,g | dJnodei) *° Gaussian where 
the observed distance has an associated uncertainty of 5 
%. 

Combining all galaxies, we have, 

prob(data | model,a,/3) = J|prob(d;bs I c^modeb /?) 

Since we do not know a or /?, we marginalise by inte- 
grating over a plausible range of both with a suitable prior, 
i.e. we integrate over 0.1 < a < 0.3 with a flat prior and 
0.1 < /3 < 0.3 with a 1//3 prior. 

Figure 2 shows likelihood contours for two different val- 
ues of Oo- 




hood analysis does not make a lot of difference to the con- 
tours whereas leaving out M31 severely limits the state- 
ment that can be made about the mass of the system. In 
the Oo = 0.99 universe, for example, we can only say that 
the combined mass of the system is probably less than 
9 X IQI^Mq. 

Figure 3 shows the same contours as figure 2 but for the 
case of an older universe. The masses generally decrease 
in this case, because less mass is needed to produce the 
same results when gravitation works over a longer period 
in time. 




M in 10"M„ 



Figure 3. As for figure 2 but with n = 1.5. 

Both figure 2 and 3 were produced from a set of solu- 
tions using 8 basis functions. We experimented with dif- 
ferent numbers of basis functions but the results remain 
the same. 

Our results indicate that most of the physics are in- 
cluded in our model if not all of them. Hence we opted for 
a likelihood analysis. The analysis shows the system domi- 
nated by M31, which makes our model not so very different 
from the one used for the timing argument. However, if 
we leave M31 out of the analysis, the resulting contours 
arc at least consistent with figures 2 and 3. There seems 
to be some slight preference for mass to cluster with the 
galaxies, but the contours are not really conclusive. 



M in 10'^ M„ 



Figure 2. Likelihood contours for variational solutions asso- 
ciated with different mass distributions: above Q. = 0.4, below 
n = 0.99 In both 1. Contour levels differ by a factor 

of 10. 

Masses are particularly high for low because the 
quadrupole force, which pulls the Milky Way and M31 

apart is higher in those cosmologies and so we require more 
mass to keep the two galaxies at the same distance from 
each other. 

The likelihood contours are dominated by M31: leav- 
ing one or more of the dwarf galaxies out of the likeli- 



4. DISCUSSION 

The A^-body check proves that the code works and is 
ready to be extended to larger systems. The main prob- 
lem with a system like the Local Group is the occurrence 
of multiple solutions. We have no guarantee that the solu- 
tions we find are the real ones, even if they fit all redshifts 
perfectly and predict distances which are not too far out 
from the measured ones. In fact, some tests suggests that 
at least for certain masses, several possible solutions are 
very close to each other so it is easy to pick the wrong 
one. This problem will not occur in larger systems, so 
they should in fact be easier to deal with. 
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In future work, two issues will need to be addressed: 
First, we need an approximate and more efficient way of 
doing the force calculations in equations |ll| and The 
direct sum that we have used needs to be replaced by a 
standard A^-body method. Second, we need more efficient 
convergence. As mentioned above, the main difficulty lies 
in preventing the coefficients from finding two-cycles in- 
stead of fixed points of the iteration. A faster yet still 
robust method for dealing with this problem is needed. 

Even when these problems are solved, the variational 
calculations will never have as many particles as A^-body 



simulations. The reason for pursuing this method is that 
unlike the A^-body calculation, which can only reproduce 
the current state of the system in a statistical sense, the 
variational method can fit the current state exactly. 
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